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Context. Eccentric accretion disks in superoutbursting cataclysmic and other binary 
^ ' systems. 
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Aims. We study the development of finite eccentricity in accretion disks in close bi- 
nary systems using a grid-based numerical scheme. We perform detailed parameter 
studies to explore the dependence on viscosity, disk aspect ratio, the inclusion of a 



> 

, mass-transfer stream and the role of the boundary conditions 

Methods. Using a two-dimensional grid-based scheme we study the instability of ac- 
cretion disks in close binary systems that causes them to attain a quasi-steady state 



with finite eccentricity. Mass ratios 0.05 < q < 0.3 appropriate to superoutbursting 
cataclysmic binary systems are considered. 

Results. Our grid-based scheme enables us to study the development of eccentric disks 
for disk aspect ratio h in the range 0.01 — 0.06 and dimensionless kinematic viscosity 
v in the range 3.3 x 10" 6 - 10~ 4 . Previous studies using particle-based methods were 
limited to the largest values for these parameters on account of their diffusive nature. 
Instability to the formation of a precessing eccentric disk that attains a quasi-steady 
state with mean eccentricity in the range 0.3 — 0.5 occurs readily. The shortest growth 
times are ~ 15 binary orbits for the largest viscosities and the instability mechanism 
is for the most part consistent with the mode-coupling mechanism associated with the 
3:1 resonance proposed by Lubow. However, the results are sensitive to the treatment 
of the inner boundary and to the incorporation of the mass-transfer stream. In the 
presence of a stream we found a critical viscosity below which the disk remains circu- 
lar. 

Conclusions. Eccentric disks readily develop in close binary systems with 0.05 < q < 0.3. 
Incorporation of a mass-transfer stream tends to impart stability for small enough vis- 
cosity (or, equivalently, mass-transfer rate through the disk) and to assist in obtaining 
a prograde precession rate that is in agreement with observations. For the larger q 
the location of the 3:1 resonance is pushed outwards towards the Roche lobe where 
higher-order mode couplings and nonlinearity occur. It is likely that three-dimensional 
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simulations that properly resolve the disk's vertical structure are required to make 
significant progress in this case. 

Key words, accretion disks - binary stars - hydrodynamics 



1. Introduction 

SU Ursae Majoris systems are a class of cataclysmic variable that in addition to undergo- 
ing ordinary outbursts undergo less frequent larger superoutbursts. During these, their 
light curves show superhumps which have a period generally slightly l onger than the 
orbital period (for a discussion of relevant observations see for example 



Patterson et al 



Ak et al 



2002 



2005). 



A natural explanation of this phenomenon is that it is associated with a rigidly pro- 
cessing eccentric accretion disk. Superoutbursts and superhumps are found to be limited 
to cataclysmic binaries with mass ratio q < 0.35. Assuming that the difference between 
the superhum p period and th e orbit al period represents the solid body precession of an 
eccentric disk. 



Patterson et al 



(|2005l ) discuss the establishment of the expected correlation 
of this period with q. 

Studies of the light curve and effects on it associated with the mass-transfer stream 
impact ing the disk have enab led an eccentricity of 0.38 ±0.1 to be found for the disk in OY 



Hessman et al 



([19921 ) , who also argue that the normal superhumps are produced by 



Car by 

tidal stresses after the maximally extended disk edge approaches the companion. Simila r 



Rolfeetal 



( 20011) 



results for IY UMa, which has a similar mass ratio, were obtained by 
These authors also noted the potential importance of the impacting accretion or mass- 
transfer stream on the eccentric disk with its varying outer radius for understanding the 
late superhump light. 

These observations motivate an interest in developing good numerical schemes for sim- 
ulating accretion disks undergoing tidal interaction with a secondary binary component. 



trie d isks using particle-based methods (e.g. 



Although, starting with I Whitehurstl (jl988), there have been many simulations of eccen 



Kunze et al 



1997 



Murray 



1998; 



Smith et al 



20071 and references therein) there has been almost no grid-based work to date that is 
applicable to close binary systems, even though such an approach, apart from constituting 
a check, is promising in the context of being able to gain access to a wid er region of param - 
eter space that may be difficult to access with particle-based methods. iHeemskerkl (|1994l ) 
performed grid-based simulations, but without including disk viscosity and was only able to 
obtain an eccentric disk when only the m = 3 compo nent of the co mpa nion's tidal potential 



was re tained. Subsequent grid-based simulations bv lStehlej (|l999| ) and lKornet fc Rozvczka 



2000) also did not show clear evidence for the onset of an eccentric instability. 



Papaloizou 



2005) solved the linearized equations for the m = 1 mode of a free disk without a binary 



companion and considered the possible development of a parametric instability when the 
vertical structure was incorporated. 
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A study of eccentric disks also has a pplication to other area s of astrophysics where disks 
and orbiting companions interact. Thus lKlev fc Dirksenl ([20061 ) have studied eccentric disks 
in the context of the very small mass ratios appropriate to disk-planet interactions where 
their existence may have important consequences for theories of planet formation. 

Motivated by these considerations, the aim of this paper is to make an extensive study of 
the dynamics of two-dimensional disks that lie in the orbital plane of a binary system, which 
consists of a primary star and a lower-mass secondary star in a circular orbit. The disk is 
centred on the primary star in a binary system with mass ratio in the range 0.05 < q < 0.3 
thought to be relevant for superoutbursting cataclysmic binaries. We perform grid-based 
simulations of disks that are initialised with circular motion, in order to study the onset of 
an eccentric instability and the properties of the resulting quasi-stationary eccentric disk. 
Use of a grid-based rather than a particle-based method enables consideration of disks 
with a wider range of viscosities and more realistic aspect ratios in the range 0.01 — 0.05 
that are presently difficult to acce ss for particle-base d methods such as SPH due to their 
inherently diffusive nature (but see lSchafer et al. Il2004l) . We establish numerical convergence 



and explore the dependence of the results on the physical parameters of the disk. We also 
consider variations of the inner boundary conditions and the presence of a mass-transfer 
stream. Due to the global nature of the eccentric disk mode, results are found to be sensitive 
to these. When a stream is present we are able to find a transition to stability when the 
viscosity is reduced below a critical value. We also find that incorporation of a mass- 
transfer stream is important for obtaining values for the disk precession frequency that are 
compatible with observations. 

The plan of the paper is as follows. In Sect. [3] we describe the hydrodynamic model. We 
then go on to describe the numerical methods we use in Sect. 12.21 Numerical simulations 
showing the growth of eccentricity for a typical disk model are then presented in Sect. [3] It 
is found that a quasi-stationary eccentric disk with mean eccentricity in the range 0.3 — 0.5 
is often produced. We explore the dependence of this phenomenon on disk viscosity, aspect 
ratio and the inner boundary condition in Sect. 01 showing that the development of an 
eccentric disk is favoured at high viscosity but, because it is a global phenomenon, it is 
sensitive to the treatment of the disk inner boundary. The quasi-stationary eccentric disk 
is found to precess as a rigid body. In this section we also study the dependence of the 
disk precession rate on the disk aspect ratio and boundary conditions showing that the 
precession rate is prograde for cool enough disks. We then investigate the dependence of 
the tendency of the disk to become eccentric on the binary mass ratio in Sect. [5] and discuss 
numerical resolution demonstrating the numerical convergence of our results in Sect. [5] 

In Sect. [7] we explore the instability mechanism and establish the importance of the 
777 = 3 component of the tidal p otential for cau sing the instability through the mode- 
coupling mechanism proposed by iLubowl (|1991af ) . Comparison with a linear analysis of 
eccentric modes is made in Sect. |5] 

We go on to investigate the effect of a mass-transfer stream in Sect. [S] This is shown to 
make an important contribution to the disk precession rate and also to cause the disk to 
become stable to developing eccentricity for sufficiently small viscosity. Finally in Sect. [TO] 
we summarize and discuss our results. 



4 Kley at al.: Eccentric disks in binary stars 

2. Hydrodynamical model 

To model a flat two-dimensional disk we generally use a cylindrical (r, tp) coordinate system 
centred on the primary star that corotates with the binary. The latter is assumed to have 
a circular orbit and thus to rotate uniformly. This coordinate system is convenient because 
the positions of the two stars remain fixed during the simulations. We have checked that 
the same results are obtained when simulations are run in the inertial frame. The disk 
material is treated as a viscous fluid with a constant kinematic shear viscosity and (except 
where noted below) vanishing bulk viscosity. 

For a flat disk located in the orbital plane z = 0, the velocity components are u = 
(u r , u v ,0). Thus u = u r is the radial velocity and the local disk angular velocity is VL = u v jr 
as measured in the corotating frame which rotates with the constant angular velocity, co, 
of the binary. 

The vertically integrated continuity equation is 

f)Y 

— + V(Su)=0. (1) 

The vertically integrated radial component of the equation of motion is 

^2 + V(£ Ur u) = Sr(, + fi) 2 -^-Ef- E6 r + f r , (2) 
at Or Or 

and the azimuthal component can be written in the form 

glgr!|±^ + V[£r> + n)u] = -| - E§£ - r£6, + r U . (3) 
Here £ denotes the surface density 
pdz, 

where p is the density, p the vertically integrated (two-dimensional) pressure. The gravi- 
tational potential $ is due the primary with mass Mi and the secondary having mass M 2 
and is given by 

GMi GM 2 
$ — 

|r-ri| |r — r 2 | 5 

where G is the gravitational constant and ri and r 2 are the position vectors of Mi and M 2) 
respectively. Since the mass of the disk in cataclysmic variables is only about 10 _10 Mq 
we can safely neglect the self-gravity of the disk as well as its influence on the binary 
orbit. The additional terms b = (b ri b v ) take into account the acceleration of the origin 
of the coordinate system. Rotation of the coordinate axes is taken account of through the 
terms involving to. The visc ous force pe r unit area is f = (/ r , f v ). The explicit form of the 



components of f is given in iKlevi (|1999f l . 

For simplicity we use a locally isothermal equation of state where the vertically inte- 
grated pressure p is related to the surface density and sound speed through 

P = Zc 2 s . (4) 

The local isothermal sound speed c s is taken to have a fixed dependence on radius and is 
given by 

£ s — — ""Kcp = hu Kcp , (5) 
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where uxcp = y GM\ / r denotes the Keplerian orbital velocity of the unperturbed disk 
that would exist in the absence of the binary companion. Equation ([5]) follows from vertical 
hydrostatic equilibrium. The constant ratio h of the vertical height H to the radial distance 
to the primary, r, is taken as a hxed input parameter. Here we adopt a standard value of 



h = H/r = 0.05 



(6) 



This is s omewhat large wh en compared to values expected for disks in cataclysmic binary 
systems (jSmith et al.ll2007j ). but smaller values are also considered below. 



2.1. Dimensionless units 

We adopt a system of units for which the unit of length is taken to be the orbital separation 
of the binary, abin- The unit of time, to, is taken to be the reciprocal of the orbital angular 
frequency oj of the binary. Thus 



1/2 



G(M 1 +M 2 ) 
The orbital period of the binary is then given by 



orb 



2vrf . 



(7) 



(8) 



The results of calculations will usually be expressed as a function of the evolution time 
in units of P or b- The unit of velocity is given by Uq = abinAo- Setting the gravitational 
constant G to be unity, Eq. ([7]) can be also viewed as specifying the unit of mass to be the 
total mass of the system M = M\ + Mi- The surface density can be scaled by any constant 
value So, because this drops out of the governing equations. Accordingly the surface density 
may be scaled to make the total disk mass be any specified value. The kinematic viscosity 
coefficient, v, is expressed in units of vq = ObinMo = wa bin- 



2.2. Numerical methods 

The governing Eqs. (H}-© determining the evolution of S, u and Q are solved using an 
Eulerian finite-difference scheme. The computational domain [r m ; n , r roax ] x [(p m in> ^ mM ] is 
subdivided into N r x N v grid cells that are equally spaced in each coordinate direction. 
For typical runs we use N r = 200, N v = 200, though many have been checked at the higher 
resolution N r = 300, N v = 300. 

The numerical method utilises a staggered mesh, where scalar quantities are located 
at cell centres and vector quantities at cell interfaces. We adopt an operator-split scheme 
for which advec tion is based on the spatially 2nd order monotonic transport algorithm of 
van Leerl (|1977l ). and the viscosity can be treated either explicitly or implicitly. The basic 
features of the two-dimensional code RH2D that we use have been described in more detail 



in 


Kiev 


( 


1989) 


Kiev 


( 


1999). 



Kiev ( 19891 ) . The use of the code in cylindrical (r, ip) coordinates has been described in 



Implementation in a rotating coordinate system re quires spec ial treatment of the 



Coriolis terms to ensure angular momentum conservation l|Klevlll998f ). To prevent possible 



instabilities due to shocks, that could for example be generated by the tidal perturbation 
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of the secondary, we use an artificial bulk viscosity as described in iKlevi (|1999f ) . This is 
in addition to the physical shear viscosity v. We note, that an excellent test of the cor- 
rect implementation of viscosity in numerical simulations of accretion dis ks is the study 
of no n-axisymmetric instabilities in the viscously spreading ring problem (jSpeith fc Kiev 



2003). 



In addition, to prevent numerical instabilities due to very low densities (or very steep 
density gradients) near the outer boundary of the disk, where it is tidally truncated by 
the secondary, we introduce a 'density floor', Sfl oor . This is applied such that whenever 
evolution causes the density at any location to fall below Sfl oor , it is reset to this value. 
Formally, this implementation does not conserve mass exactly, but for sufficiently low values 
it does not produce any influence on the dynamics. In all simulations we use a value of 
Sfloor = 10~ 8 expressed in units such that the initial value of £ at r — a D in would be unity 
in the absence of tidal truncation (i.e. in units of Eo defined in Eq. [TU] below). We test the 
dynamical influence of the floor in section 14.31 below. 

Since this is an explicit code, the timestep limitation is given by the Courant condition 

At = / c min^M (9 ) 

grid | It | +C S 

where the Courant number f c has to be smaller than one, and we typically use f c — 0.70. 
After a few test simulations it soon became clear that, due to the large range in radius 
with r max /r m ; n typically larger than 10, the required computational times become excessive 
because of the stringent timestep limitation arising from the inner boundary. This occurs 
because the Keplerian angular velocity increases rapidly as the primary star is approached 
according to Qk oc r~ 3 / 2 . This results in approximately 25,000 timesteps being required 
for one binary orbit with our standard resolution. As typically a few hundred orbital 
times need to be calculated per run, this necessitates about 10 million timesteps for each 
run. Exploration of a large para meter space w ould be prohibitive. For that purpose we 
implemented the FARGO algorithm (|Massetll2000l ) into our RH2D code. In this algorithm, for 
each annular grid ring, the mean angular velocity Cl is calculated and all variables are simply 
shifted by an integer number of grid cells, corresponding to an advection velocity as close 
as possible to rCt. Advection is then only performed using the residual velocity, leading to 
a substantially increased timestep. This method leads, for our standard parameters using 
an equidistant grid, to a speed-up by a factor of about 5. 

Additional test simulations have shown that for cases that do not use the FARGO algo- 
rithm, numerical instabilities due to very small densities are likely to occur in the outer 
region. These could only be avoided by further reducing the timestep or adopting an im- 
plicit treatment of viscosity, both of which result in additional computational cost. A higher 
density floor also helps to reduce these instabilities but to avoid them entirely (and still 
use a large timestep) one would need unreasonably large values that possibly influence the 
dynamics; see the discussion on boundary conditions in Sect. POl below. 

Here, the use of the faster advection algorithm enhances the stability properties, allow- 
ing an explicit treatment of viscosity. A few comparison simulations that do not use FARGO, 
but allowed a redu ced timestep, were performed to check the accuracy of our method using 
the NIRVANA code (|Zieglerlll998l ). Finally, we note that the order of the applied numerical 
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scheme matters in obtaining the eccentric state of the disk. Switching to a simple first order 
upwind scheme does not result in eccentric disks no matter how fine the grid resolution is 
chosen. For the 2nd order upwind scheme the different choices of the slope-limiter all yield 
very similar results. 




X 

Fig. 1. The grid structure used in the computations. The radial domain is [r m i n , r max ] — 
[0.05,0.70] and the tp domain is the whole annulus [0, 2tt]. For our standard case this is 
covered with N r x N v = 200 x 200 grid cells that are equally spaced in each direction. The 
Roche lobe of the primary, for our standard mass ratio q = M2/M1 — 0.1, is defined by 
the additional solid (red) line. The primary and secondary star are located at (x — x p — 
0,2/ = Up — 0) and (x — x s = —l,y = y s = 0) respectively. Note that in this figure, for 
illustrative purposes, we display a low-resolution 50 x 50 grid. 



2.3. Setup and boundary conditions 

The coordinate system is centred on the primary star and, because it corotates with the 
binary system, the position of the secondary in Cartesian coordinates is fixed at x = 
— l,y — (see Fig. [1]). For the radial domain we adopt [r m j n , r max ] where r m i n — 0.05 and 
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r max = 0.70 for all models in units where the binary separation is unity. The azimuthal 
domain is taken to be [0, 2tt]. 

For the initial surface density profile, for r < 0.4 in dimcnsionless units, we adopt a 
power law of the form 

E(r,p) = £ r- 1/2 . (10) 

For r > 0.4, £ is reduced smoothly towards a minimum value of 10 _8 So- The initial velocity 
is taken to be pure Keplerian rotation around the primary star as would be measured in 
an inertial frame (u r = 0, ril — {GMi/r) 1 / 2 — ru>). No correction is made for pressure 
effects or the non-axisymmetric binary potential. The presence of the secondary star is 
allowed to cause the surface density distribution to evolve away from the initial profile 
on a dynamical timescale. This normally causes a rapid adjustment of the outer profile 
through tidal truncation of the disk. The fixed temperature stratification is always given 
by T(r) oc r _1 which follows from the assumed constant aspect ratio h = H/r = constant. 

In the azimuthal (ip) direction, periodic boundary conditions for all variables are im- 
posed. At the outer radial boundary of the grid which extends somewhat beyond the Roche 
lobe of the primary (see Fig.[T]) we use standard outflow boundary conditions. This imposes 
zero gradient conditions on all variables apart from the radial velocity, i.e. the values in the 
ambient ghost cells have the same values as the last active grid point. The radial velocity 
at the ambient ghost cell is set equal to the value at r — r max if the flow direction is out- 
ward (zero gradient condition) and zero otherwise. This outflow condition allows material 
to leave the system but inflow is not permitted. 

As we found that the form of the boundary condition applied at the inner boundary 
may play an important role in determining the characteristics of the flow, we implemented 
four different boundary conditions at r m j n : 

— Rigid reflecting: Here the inner boundary acts as a solid wall where the radial velocity 
is set to zero, while for the other variables mirror symmetry conditions are applied to 
determine values in ghost cells. 

— Viscous outflow. Here an axisymmetric form for the radial velocity at the inner boundary 
is specified. This is directed inward with a magnitude given by 

3 is 

u(r min ) = --- . (11) 



nun 



This boundary condition imposes a steady-state accretion disk flow. 

Damping condition: All variables within a certain region near the inner boundary are 

relaxed on the orbital timescale at r m \ n towards their initial val ues. This approach 



de Val-Borro et al 



reduce s wave reflection from the inner boundary and is described in 
(2006). In addition u r (r m i n ) is fixed at zero for all time. 
— Free outflow: This boundary condition corresponds to the outflow condition applied at 
'"max- The radial flow velocity at r m i n can only be directed inwards and no inflow from 
the inner regions into the computational domain is allowed. 

We found from our numerical simulations that the form of the inner boundary condition 
can make quite a difference in the results concerning the eccentricity and precession rate 
of the disk. 
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3. Standard model 



9 



To study the effect of varying the physical and numerical parameters defining the simula- 
tions, we first consider in detail a standard model for use as a reference for the extensive 
parameter studies we have performed. The parameters defining this model are summarized 
in Table [TJ 



Mass ratio (q = M 1 /M 2 ) 


0.1 


Disk thickness (h = H/r) 


0.05 


Viscosity {v) 


10~ 5 


Radial range ([r min , r max ]) 


[0.05,0.70] 


Grid {N r x N v ) 


200 x 200 


Inner boundary 


reflecting 



Table 1. Parameters of the standard model. Values are quoted in dimensionless units. 




Fig. 2. Greyscale plot of the surface density of the disk for the standard model after 1 
binary orbit. The solid (red) curve indicates the Roche lobe of the primary star (central 
red dot). 



10 
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Fig. 3. Greyscale plot of the surface density of the disk for the standard model at 
100, 150, 175 and 200 orbital periods. The solid (red) curve indicates the Roche lobe of 
the primary star (central red dot). 

We now describe the characteristic behaviour of the standard model. Due to the strong 
gravitational perturbation induced by the secondary, the initial axisymmetry of the disk is 
rapidly destroyed and a strong two-armed spiral density response is formed. This is shown 
in Fig. [3 where the surface density is displayed after the system has evolved for one orbital 
period of the binary. The two-armed spiral which is strongest in the outer disk region but 
extends all the way through the disk down to the primary star is clearly visible. The surface 
density shows some evidence for a non-axisymmetric pattern near the inner boundary of 
the disk, confined to the innermost two radial gridcells. Such a feature (here a m = 6 mode) 
is generated typically in grid-based simulations with reflecting inner boundary conditions, 
and is not a consequence of using the FARGO-algorithm. It can be reduced or minimised 
by using special, non-reflecting or damping, boundary conditions. After the onset of the 
eccentric instability this feature dissapears in the simulation. 

In the subsequent evolution the disk slowly becomes eccentric. This is indicated in 
Fig. [3] where we display the two-dimensional surface density structure at four different 
times. The distortion of the disk is evident at later times. While in the first snapshot taken 
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after 100 orbits there is only a slight indication of disk eccentricity, the subsequent panels 
demonstrate a growing eccentric disk mode and one can infer the presence of an apsidal 
line joining the locations of minimum and maximum displacements from the primary. 
Additionally there is indication for precession of this apsidal line. 



10 
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Fig. 4. Radial distribution of the azimuthally averaged surface density £(r) at different 
times. The solid black curve indicates the tapered initial profile. 

The tidal torques induced by the secondary also lead to a truncation of the disk and a 
rearrangement of the (mean) radial surface density profile. This is demonstrated in Fig. [4] 
where the azimuthally averaged density is shown at 5 different times. The solid black 
curve indicates the initial (axisymmetric) profile, £ = r -1 / 2 , with a taper applied for 
r > 0.4. Clearly, the disk becomes truncated at a radius around r « 0.5 where the surface 
density starts to drop much more steeply. After about 250 orbits the average surface density 
distribution attains a new quasi-stationary profile (see below). 

In Fig. [5] we display the time-evolution of several global disk quantities (in dimensionless 
units). The total mass in the computational domain is seen to drop off rapidly within the 
first 2 orbits due to the initial mass loss through the outer boundary induced by the tidal 
perturbation of the secondary (see also Fig. [5]). After this initial burst of mass loss, the disk 
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Fig. 5. Evolution of various global disk properties of the standard model as functions of 
time. The total mass in the computational domain, or equivalently the disk mass, is shown 
in the upper left panel. The disk radius, defined as the radius containing 99% of the mass, 
is shown in the upper right panel. The disk eccentricity is shown in the lower left panel 
and the m = 1 and m = 2 Fourier amplitudes are shown in the lower right panel. 



is tidally truncated, such that as seen in Fig. the disk radius, defined to be the radius 
containing 99% of the mass, is reduced from 0.61 to 0.48. This is then smaller than the 
Roche lobe. The mass subsequently remains constant until t w 200. The disk radius reaches 
a minimum at around t 80, and increases slowly thereafter. This increase of the radius of 
the disk is coupled with an increase in the mean eccentricity ed of the disk which is displayed 
in the bottom left panel of Fig. The quantity ed is determined by first calculating the 
eccentricity for each individual gridcell assuming that it is in an unperturbed orbit around 
the central star, with the position and velocity given by the actual grid values averaged to 
the cell centres. To obtain the global value ed we then perform a mass-weighted average 
over the individual cell values. The longitude of the disk's p eriastron zu^ is calculated in 
the same way. The explicit formulae are given for example in iMurrav fc Dermottl ([1999) . 

As seen in Fig. [SJ ed begins to exponentiate after t s=s 20 with an e-folding time of about 
80 orbits. After about 250 orbits the eccentricity has settled to a final value of ed = 0.37 
with an average outer disk radius rd ~ 0.55. The increase in ed is also associated with a 
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further decrease in the disk mass. Beyond t w 250 mass spills over the Roche lobe and is 
lost from the system. This mass loss does not affect the dynamics of the system (i.e. it does 
not change the equilibrium values of ed and rj) because the absolute value of the mass 
density scales out of the evolution equations. This feature has been confirmed numerically 
in additional test calculations where the total mass in the system is kept constant by 
rescaling the density of all grid cells at each timestep in order to account for the lost mass. 
In the bottom right panel of Fig. [5] the time-evolution of the m = 1 and m = 2 surface 
density Fourier amplitudes calculated at a single disk radius (r = 0.37) are displayed. The 
m = 1 mode associated with the eccentricity increases exponentially from the beginning 
of the run up to t ss 200 after which induced additional mass loss leads to a reduction in 
the amplitude. To show the evolution more clearly we have plotted here only the first 200 
orbits. The only other Fourier amplitude that plays a significant role is that associated 
with m — 2, which at early times (when the spiral arm feature dominates the disk) is much 
larger than the m — 1 amplitude, but it becomes weaker at later times. In the final stages 
it increases again, saturating eventually at a value still smaller than the m = 1 amplitude. 



As mentioned above, Fig. [3] indicates global precession of the eccentric disk. This is 
quantified in Fig. [5] where the solid (red) curve corresponds to the time-evolution of the 
longitude of pericentre (periastron) of the disk, Wd, calculated in the same manner as e<j- 
Obviously the disk precesses as a whole at a constant rate, a feature that is supported by 
the additional light (green) curve which corresponds to the phase of the m — 1 Fourier 
component at one given radius r = 0.37. Both angles precess at the same rate, confirming 
our interpretation. The periodic disturbances in ra7<j coincide with the times when the disk's 
apocentre is closest to the secondary star. In this case the outer parts of the eccentric disk 
are strongly disturbed which renders the calculation of Wd more inaccurate. From the plot, 
which shows the phases in the inertial frame, we conclude that the disk experiences a slow 
retrograde precession with a period « 16.4 orbits. 



To study the internal structure of the disk we show in Fig.[7]the disk's radial distribution 
of eccentricity and the longitude of pericentre at different times (indicated by the legend) 
which lie within the time interval plotted in Fig. [5] The distribution of eccentricity within 
the disk does not vary much with the binary orbital phase and drops smoothly towards zero 
in the centre of the disk, due to the rigid boundary at the inner radius. In the main part 
of the disk the eccentricity lies around e = 0.4 which is consistent with two-dimensional 
surface plots of the density (see the snapshot at t = 200 in Fig. [3]). The radial variation 
of the pericentre indicates a small twist within the innermost parts of the disk, which is 
larger when the disk's apocentre is closest to the secondary star (at times 351 and 352) 
when it can reach about 45 degrees. In addition, the fact that the curves at different times 
appear vertically shifted indicates again the precession of the disk as a single entity. 
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Fig. 6. Time-evolution of the longitude of pericentre (periastron) of the disk. The solid 
darker (red) curve is the mass-weighted mean pericentre over the whole disk, calculated 
as described in the text. The lighter dashed (green) curve corresponds to a phase of the 
m = 1 Fourier mode at the single radius r = 0.37. Note that the phases are calculated in 
the inertial frame. 

4. Dependence on viscosity, temperature and the inner boundary condition 

4,1. Viscosity 

The first parameter we investigate is the kinematic viscosity v. To study dependence of the 
results on the magnitude of the kinematic viscosity we started from our standard model 
(y = 10 ) and varied only the value of v keeping all other parameters unchanged. 

The variation of the mean disk eccentricity with time is displayed for different kinematic 
viscosities in Fig. [51 where the middle (green) curve refers to the standard model which 
(for reference) has been displayed already in Fig. [5] at early times. The growth rate of the 
eccentricity depends strongly on the magnitude of the viscosity: increasing the viscosity 
shortens the growth time substantially, while for smaller viscosities it takes much longer 
for the disk to become eccentric. For the lowest- viscosity model (v = 3 x 10~ 9 ), which 
corresponds essentially to an inviscid calculation, the disk does not attain significant ec- 
centricity in a run time of 1000 orbits. However, there is an indication that there may be 
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Fig. 7. Radial variation of the eccentricity and longitude of pericentre of the disk at various 
times (quoted in binary orbits). 
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Fig. 8. Time-evolution of the mean eccentricity of the disk (natural logarithm) for different 
kinematic viscosities (in dimensionless units). 

a very long-term eccentricity growth produced by numerical effects. From the upper left 
panel of Fig. [5] a simple extrapolation indicates a limiting value of v for which the disk be- 
comes eccentric of v ~ 2 x 10~ 6 . Already the growth time for v = 3.3 x 10~ 6 is ~ 162 orbits. 
In general the models settle to a quasi-stationary state with a constant eccentricity, which 
decreases as the viscosity decreases. However, quasi-steady states with mean eccentricities 
below e = 0.36 do not seem to be reached for this model by decreasing the viscosity alone. 
But then a longer time is required for this value to be attained. 
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Fig. 9. Dependence of various global disk parameters on the magnitude of the viscosity. 
Upper left panel: Growth rate of the eccentricity of the disk; Upper right panel: The radius, 
rd, within which 99% of the mass of the disk is contained; Lower left panel: The mass- 
weighted mean eccentricity of the disk ed; Lower right panel: Precession rate of the disk 
Wd, negative values indicating retrograde precession. The values for rd, ea and Wd are 
determined at a time at which the system reaches its final quasi-steady state. 

Clearly the initial growth of ed in these models is exponential and growth rates, ad, 
defined through 

ed(t) oc exp(tTdi), (12) 

can be determined by fitting simulation data to the above relation. The results are plotted 
in the upper left panel of Fig. [5] where the variation of ad with viscosity is plotted. In the 
remaining 3 panels we plot the equilibrium values of the radius rd , the eccentricity ed and 
the precession rate vjd of the disk. Each of the quantities (cd, ?"d, ed) shows a strikingly 
similar dependence on v. This behaviour can be understood in terms of the diffusive in- 
fluence of viscosity. For larger viscosity the disk material will tend to spread outwards at 
a faster rate with more material being pushed to the outer regions of the disk where it 
feels a stronger (eccentricity enhancing) tidal perturbation due to the secondary. Hence 
the increase of the growth rate with viscosity. For the very same reason, the equilibrium 
value of the radius rd as well as the final eccentricity will increase with v. The fact that 
all these quantities respond similarly to variations of v indicates that they respond in a 



Kley at al.: Eccentric disks in binary stars 17 




Fig. 10. Four snapshots of the disk surface density contours during one orbital period as 
viewed in the non-rotating frame. Time increases from top to bottom and left to right The 
solid (red) curve indicates the Roche lobe of the primary star. The viscosity is 10 -4 , the 
mass ratio is 0.1 and H/r = 0.05. As the companion has a close approach to the disk outer 
edge a tidal tail is pulled out (top right). This subsequently rejoins the disk (bottom left) 
which remains relatively unperturbed until the next close approach. 

similar manner to the outward expansion of the disk. On the contrary, the disk precession 
rate w& depends only weakly on the viscosity (for v Z 5 x 10~ 6 ) . From the lower two 
panels (ed, t^d) there is an indication that for v < 10~ J the disk shrinks rapidly causing 
corresponding changes to the precession rate while the eccentricity levels off. However, we 
did not consider simulations with v < 3.3 x 10~ 6 as it is likely that at this stage they 
become resolution-limited such that the numerical viscosity cannot be neglected. 

Finally we illustrate the behaviour of the eccentric disk once it has attained its final 
quasi-steady state for the case with v = 10~ 4 in Fig. [Til] Cases with other values of v 
show similar behaviour. In Fig.[TU]we show the surface density contours during one orbital 
period as viewed in the non-rotating frame. In this frame the secondary (moving counter 
clock wise) has a close approach to the slowly processing disk's apocentre approximately 
once every orbit. As the companion approaches the disk's outer edge a tidal tail is pulled 
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out. This subsequently rejoins the disk which remains relatively unperturbed until the 
next close approach. As even at apocentre, the disk matter rotates faster than the binary 
so the tidal tail leads the secondary and is clearly responsible for angular momentum 
transfer to the orbit. This must also be associated with enhanced tidal dissipation. As the 
disk appears relatively unperturbed away from these closest approaches, most of the tidal 
angular momentum transfer may occur as a result of these close approaches. 
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Fig. 11. Mean disk eccentricity ed and disk radius as functions of time for the standard 
model (red curve), the same case but with an added explicit bulk viscosity I'buik = 2^ 
(green curve) and a case with standard parameters, but for which the initial disk radius is 
taken to be rd(0) = 0.3 (blue curve). 



4.1.1. Effect of bulk viscosity 

The imposed kinematic viscosity, v 1 is intended to represent the effect s of some form of 
disk t urbulence, probably arising from the magnetorotational instability IjBalbus fc Hawlev 



1998). It is possible that an effective bulk viscosity may also result which could also act as a 



stabilizing mechanism. In addition, in typical SPH simulation s on this problem a n explicit 



bulk viscosity is included; see for example the recent study bv lSmith et alj (|2007l ). In order 
to investigate its effects, we performed runs for which an explicit bulk viscosity C, = EtWk 
was added. In this study we used a constant t-wik and present in Fig. QT] results for the 
case when the bulk viscosity coefficient is twice as large as the shear viscosity coefficient: 
^buik = 2ia The influence of the added bulk viscosity on the simulations is very small: the 
only effect is a slightly increased damping, as indicated for example by the longer growth 
time. Otherwise an added bulk viscosity has no dynamical influence. This indicates that the 
direct effect of viscosity on the eccentric mode is unimportant and implies that mechanisms 
such as viscous overstability cannot be responsible for the eccentric instability seen here. 
We also plot the evolution for a case where the initial disk radius ?"d(0) = 0.3, rather than 
the usual rd(0) = 0.5, in Fig. QT] As expected the initial growth is slowed down due to 
the small outer disk radius which remains around = 0.47 until t = 150. Only when the 
eccentric instability, which appears to grow (though more slowly) even for this small disk, 
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has developed sufficiently does the radius of the disk increase again, (ft should be noted 
that, owing to the way in which it is defined, an increase of may be associated purely 
with a growth of the eccentricity of the disk rather than an increase of the semimajor axis 
at which the disk is tidally truncated.) At later times the behaviour is in fact identical to 
the standard model. Hence, different initial conditions lead to the same outcome, as one 
might expect for viscous flows. 

4.2. Temperature 

In this section we study the effect of variations of the disk temperature or, equivalently, 
the relative scaleheight h = H/r of the disk. Starting again from the standard model we 
vary h from 0.01 to 0.06. The variation of e<i(i) is displayed in Fig. [12l where the thick 
(coloured) lines refer to different values of h, and the dashed lines are eyeball fits made 
during the initial growth period. The growth timescale depends on the disk temperature 
and it is found to be shorter for cooler disks. The final disk eccentricities do not depend 
on the value of h as long as h > 0.03. For the smallest value displayed here, h = 0.02, the 
standard grid resolution (200 x 200) is not sufficient. In additional simulations we found 
that runs (for h — 0.02) with higher resolution yield larger values of ed more consistent 
with the other runs. The final radius of the disk, r^, also does not depend on the disk aspect 
ratio for h > 0.02. All runs yield an outer disk radius « 0.55, with a slight tendency for 
larger disks to be associated with larger h. 
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Fig. 12. Time-evolution of the mean eccentricity of the disk for different disk aspect ratios 
h = H/r. 
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Fig. 13. Left: Dependence of the growth time (crj 1 ) of disk eccentricity on h = H/r. 
Quoted is the time (in orbital periods) to increase the eccentricity by a factor of e; Right: 
The equilibrium surface density profile of the disk normalized to the central value for 
different h. 



The inferred growth times as a function of h — H/r are plotted in the left panel of 
Fig. 021 The growth is fastest for cooler disks and the growth time l/oy increases linearly 
with increasing h. For the thinnest disks we study with h = 0.02, the e-folding time is 
around 10 orbits (for the standard viscosity). In the right panel of Fig. [TS] the equilibrium 
surface densities are displayed for various h. Obviously the value of h does not play a role 
in determining the equilibrium structure. At higher resolution, the surface density profile 
of the cool h — 0.02 disk is in better agreement with the others as well. 

We plot the equilibrium precession rate of the disk as a function of h in Fig. [I4j The 
inverse of this quantity gives the period of one complete turn of the disk (in the incrtial 
frame) . Values obtained by solving the linear free m = 1 mode eigenvalue problem using the 
mean surface density distributions shown in Fig. ll3l an d assuming a fr ee outer boundary are 



Goodchild fc Qgilvie 



also p lotted (for a discussion of this approach see lPapaloizoul (|2005l ) or 
(|2006l )). These are in good agreement with the results derived from the nonlinear simula- 
tions. For very cool disks, h £> 0.025, the precession rate is positive while for hotter disks 
the disk precesses in a retrograde sense. The precession period at the transition is infinite, 
corresponding to the disk being stationary in the inertial frame. In general the precession 
rate is determined by a combination of different effects. There is a contribution, n7d yn , 
due to the axisymmetric component of the binary potential, and there is also a pressure 
contribution, ri7 pross . The first part Wd yn is positive, giving a prograde contribution, while 
the latter part ri7 pross is negative, tending to produce retrograde precession. Our results 
confirm the existence of these contributions and indicate a parabolic dependence of Wd 
on h for locally isothermal disks. The influence of the viscosity on the precession rate is 
negligible as we have indicated above. 

A change of the direction of precession sometimes occurs during the earlier growth 
phases of the simulations and is illustrated in Fig. [T5] where we plot the time-evolution of 
the apsidal line of the disk for h = 0.03. Initially, during the period of eccentricity growth, 
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Fig. 14. Dependence of the equilibrium disk precession rate rod on disk temperature. The 
two separate lines with the triangles refer to two different grid resolutions, the upper to 
the standard case and the lower to 300 x 300. The additional curves with the small open 
and solid diamonds are approximate parabolic fits to the data points. Values obtained 
by solving the linear eigenvalue problem using the azimuthally averaged surface density 
profiles are also indicated by the filled circles. The open triangles are obtained from models 
where the disk is impacted by a mass-transfer stream (see Sect. [9] below). 

the disk precesses in a prograde sense and at time t ~ 130, as it settles into a quasi- 
stationary state, it turns around and undergoes retrograde precession. Thus this value of 
h lies close to the borderline separating these two possibilities for the final state. 

We also remark that the precession rate is sensitive to the mass distribution in the 
outer regions of the disk and it can be changed from retrograde to prograde by supplying 
matter from the secondary through the L\ point (see Fig. [T4l and Sect. [3] below). 



4.3. Effect of the inner boundary condition 

We here discuss the effects of the different kinds of boundary conditions that we have 
implemented as described in Sect. 12.31 and discuss briefly the influence of a higher density 
floor. The evolution of the mean eccentricity of the disk for the different inner boundary 
conditions for a kinematic viscosity of v — 10~ 4 is shown in Fig. 1161 

The results shown in Fig. \W\ indicate that the adoption of an inner boundary condi- 
tion with damping properties (such as the damping boundary conditions and the outflow 
condition) lead to configurations that differ from the standard case. Both the growth rate 
and the mean equilibrium eccentricity of the disk are strongly affected. Increased damping 
yields smaller growth rates and smaller final ed- The outer disk radius shows similar ten- 
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Fig. 15. Time-evolution of the pericentre (in the inertial frame) of the disk, Wd(t), for a 
model with an intermediate H/r = 0.03. 

dency but it varies within a much smaller range from rd ~ 0.60 for the damping boundary 
condition to rd « 0.64 for the viscous outflow condition (see Tabled]). The viscous outflow 
condition and the standard reflecting condition behave very similarly because both have 
strong reflective properties. The most pronounced influence of the boundary is on the pre- 
cession rate of the disk, as is indicated in Table [H For the free outflow condition we find a 
very slow positive, i.e. prograde, precession of the disk. However, as is indicated in Fig. 1171 a 
very pronounced eccentric inner hole forms in the centre of the disk which does not happen 
in cases where the inner boundary is forced to be axisymmetric. When an open boundary 
condition is used, material on an eccentric orbit can leave the system freely and is not 're- 
flected back' into the computational domain. A side-effect of the use of the open boundary 
condition (and also the damping condition) is a reduction in the overall eccentricity of the 
disk. Clearly, the treatment of the inner boundary condition is an important issue that 
needs to be considered in more detail. In the case of a central compact stellar object, such 
as a white dwarf or neutron star, which has a very small pressure scaleheight in its upper 
layers, it might behave approximately as a rigid wall with typical reflecting wall conditions 
for the normal velocity and a no-slip boundary condition for the tangential velocity. This is 
the situation of the boundary layer problem in the theory of accretion disks which has an 
associated radial lengthscale that is too small to make inclusion tractable in a simulation 
of the type considered here. However, it appears that the overall disk properties may be 
significantly affected by the behaviour of this small innermost region of the disk. The final 
line in Tab. prefers to a model with a reflecting inner boundary and a high density floor of 
10~ 3 (instead of 10 -8 ). While the final values of the disk's eccentricity and radius are not 
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Fig. 16. Time-evolution of the mean eccentricity of the disk for different inner boundary 
conditions for a kinematic viscosity of v = 10~ 4 (ten times the standard value). 

changed much, there is a significant difference in the precession rate. The higher density 
in the outer parts of the primary's Roche lobe slows down the precession by a factor of 
about 2.8. Hence, to accurately determine the important rate of disk precession a small 
floor value is clearly required. 



Boundary condition 


ed 


T& 


cbd[l/-Porb] 


Reflecting 


0.50 


0.640 


-0.064 


Viscous outflow 


0.56 


0.651 


-0.060 


Damping 


0.31 


0.598 


-0.046 


Free Outflow 


0.33 


0.597 


+0.0037 


Reflecting, high floor 


0.52 


0.620 


-0.022 



Table 2. Results obtained using the different inner boundary conditions and density floor. 



5. Mass ratio 

We now investigate the influence of the mass ratio q — M%/Mi on the disk dynamics. To 
speed up the simulations all of the runs in this section have been performed with v = 10~ 4 
which is ten times larger than adopted for the standard simulation (see Table 1). 

The eccentricity as a function of time is shown in Fig. [T5] for different q with over- 
plotted straight-line fits indicating exponential growth. For larger mass ratios q the growth 
is strongly enhanced, but all models settle to approximately the same final eccentricity 
ed- For the smallest mass ratio {q — 0.05) there is no direct indication of an exponential 




Fig. 17. Two-dimensional surface density structure for the model with the open outflow 
inner boundary condition after 250 orbital periods. 



growth phase and we have drawn an approximate straight-line fit applicable to the middle 
of the growth phase. The extracted growth rates are plotted as a function of mass ratio in 
the left panel of Fig. [THl For small q the growth rate increases linearly with q but it slows 
down at the largest value q = 0.3 plotted. We have performed explorative runs for larger q 
and find eccentricity growth all the way to q — 1 with fastest grow occuring at q ss 0.3 (see 
Fig. [TO]). Below it is argued that this behaviour may be related to higher order non-linear 
mode coupling which becomes more important for larger q, but in this work we do not 
elaborate on this issue any further. 

The variation of the disk radius with time (right panel of Fig. shows a sharp drop 
in rd just after the start of the simulations due to the tidal effect of the secondary which is 
stronger for larger q. After attaining a minimum value, as low as 0.46 for the largest q, the 
disk expands again on the eccentricity growth timescale. The quasi-stationary disk radius 
is smaller for larger q due to the enhanced tidal effects. Note that for the q = 0.05 case we 
have used a larger outer boundary radius r max = 0.77. 
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Fig. 18. Time-evolution of the mean eccentricity of the disk (logarithmic) for different 
mass ratios q = M 2 /Mi of the binary for a kinematic viscosity of v = 10 -4 (ten times the 
standard value). The thinner black dashed lines indicate approximate linear fits. 
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Fig. 19. Left: Dependence of the growth rate (rjd) of the disk eccentricity on the mass 
ratio q. Right: The disk radius as a function of time for different q. The curves are 
smoothed. 

6. Numerical resolution 

An important issue in (multi-dimensional) simulations is the question of numerical reso- 
lution. To study its effect we have varied the grid resolution in our standard model from 
100 x 100 up to 400 x 400. The grids we use are always uniform. The models are typi- 
cally started at t = to study the whole growth phase of the evolution. In some cases 
higher-resolution simulations have been started using results obtained from coarser grids 
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Fig. 20. Left: Dependence of the growth rate (<7d) of the disk eccentricity on the numerical 
resolution. The two-dimensional grid always has an equal number of grid points in each 
direction. Thus in the standard case 200 corresponds to 200 x 200. Right: The precession 
rate of the disk as a function of resolution. 

to provide initial data with the object of investigating the quasi-stationary state. Our main 
results are contained in Fig. [20l where the dependence of the growth and precession rates 
of the disk are indicated for different grid resolutions. The physical parameters used are 
those of the standard model. At low resolution, the growth rate increases approximately as 
a linear function of the number of grid cells. Beyond 300 x 300 it levels off which indicates 
that we have reached numerical convergence. This finding is corroborated by the fact that 
the precession rate apparently levels off at the same resolution (right panel). The converged 
growth rate Cd ~ 0.016 differs from the value for our standard resolution (cm ~ 0.012) by 
approximately 30%; similarly the value of Wd differs by about 25%. While the absolute 
values for the growth and precession rate change, the qualitative behaviour of the disk, in 
particular the prograde slow precession in this case, does not change. 



7. Instability mechanism 

7.1. Mode coupling 

In this section we investigate the mechanism of excitation of the disk eccentricity. A mech- 
anism responsible for this that involves the coupling of velocity (or surface density) pertur- 
bations with different azimuth a l mode numbers through the compani on's tidal poten tial 
has been described by iLubowl (|1991a 



This was later confirmed by ILubowl (|1991a ) by 
carrying out and analysing appropriate SPH simulations. 

Although incorporating the slow precession of the disk as viewed in the inertial frame in 
the discussion below would introduce minor adjustments, for simplicity we neglect it. Then 
when viewed from this frame, the eccentricity of the disk may be considered to be associated 
with a stationary surface density perturbation with azimuthal mode number m — 1. A disk 
in which the gas orbits in circles about the primary is subject to a tidal potential with 
components which vary as trigonometric functions of m(ip — cut), m = 0, 1, 2, 3, . . . These 
produce corresponding responses in the state variables that scale with the mass ratio, q. 
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When the disk becomes eccentric through the excitation of an m — 1 mode that is 
approximately stationary in the inertial frame, nonlinear coupling between the response 
to the tidal potential and this mode produces an effective forcing which will be a sum of 
components that are trigonometric functions of (m ± l)ip — rmot, m = 0, 1, 2, 3 . . . These 
will have amplitudes that scale as eq, e being the characteristic disk eccentricity. More 
generally a coupling of order j produces effective forcing (m ± j)ip — mcot, m = 0, 1, 2, 3 . . . 
with pattern speed muj/(m ± j) and amplitude scaling as e J . 



7.2. Energy and angular momentum exchange with the companion 



Each of these effective forcings produces a surface density response that results in energy 
and angular momentum exchange with the companion. This exchange may in some cases 
produce a feedback causing the eccentricity to grow. This works because a companion in 
fixed circular orbit must undergo angular momentum changes, A J, and energy changes, AE, 
such that AJ = AE/ui. However, the angular momentum and energy changes associated 
with a disturbance with pattern speed mw/(m±j) are such that AJ = AE(mij)/(muj). 
Thus for surface density disturbances with the — alternative, that rotate faster than the 
companion and are thus expected to transfer orbital angular momentum to it, too little 
is transferred per unit of energy transferred. The difference has to come from the disk 
orbits which accordingly become eccentric. Energy and angular momentum transfer to the 
companion i s expected to be especially ef fective if there is an associated Lindblad resonance 



in the disk ( Goldreich fc Tremaine 



1980) 



Similarly disturbances with the + alternative, which rotate more slowly than the com- 
panion and are expected to gain angular momentum, may cause eccentricity growth. 
However, such disturbances are less closely match ed to the rotation of the dis k and in 
particular have no Lindblad resonances within it ([Goldreich fc Tremaine 1979 ) and are 
thus not expected to be significant. 

Accordingly, disturbances rotating with pattern speeds muj / (m — j), with m and j ^ 
m being integers, may be associated with eccentricity growth if there is an associated 
significant energy and angular momentum transfer. The presence of a Lindblad resonance 
is expected to enhance this. Although in principle any value of j might contr ibute to 



1991a ) and the 



eccentricity growth, linear instability could result only when j = 1 (jLubowl r 
forcing as a result of coupling is linear in e. In that case the only value of m leading to a 
disturbance with a Lindblad resonance within the disk is m — 3. From the above discussion 
this disturbance has azimuthal mode number m — j = k = 2 and pattern speed 3w/2. The 
associated Lindblad resonance corresponds to the 3 : 1 resonance at which the disk orbital 
period is one third of that of the binary. Wh en this is present the mode coupling may lead 
to a linear instability to eccentricity growth (jGoodchild fc Ogilviel l2006). 

According to the above view of the mode coupling leading to eccentricity growth, the 
surface density perturbation corresponding to the disk eccentricity should grow together 
with a surface density perturbation with k = 2 and angular frequency lu>, with 1 = 3. 
This disturbance is derived from a coupling between the disk eccentricity and the m = 3 
component of the tidal potential, which accordingly plays an important part in the process. 
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7.3. Mode amplitudes 



Following lLubowl (|1991bf ). in order to demonstrate the operation of this mechanism and the 
significance of the m — 3 component of the tidal potential, we define the mode amplitude 
associated with the pair (k, I) as 



Ski = VKi 



In this expression 

rt+2-rr 



o-ki = It ^ I a ^( r ' L P 1 cos (k(p) cos(luit') dA dt' 
b k i = J f t+2?r J A S(r, if, t') sin(fc^) sm{lujf) dAdt' 



cu = Jt +27r Ja ^( r ' fi cos (kf) sin(lu>t') dA dt' 



du — Jt +2 ^ J A T,(r,(p,t')sin(k(p)cos(lujt')dAdt' , 



(13) 



(14) 
(15) 
(16) 
(17) 



where the integration is over the surface area of the disk, A. Note that Su also involves 
integration over the domain (t, t + 2?t) and is thus a function of the local time t. 

In order to test the predictions of the above mode-coupling analysis, 
we evaluate the time-evolution of the mode amplitudes for the (k, I) pairs 
(1,0), (1,3), (2, 3), (3, 3), (2, 2), (3, 2) for runs with q = 0.1 and q = 0.2. We describe 
below results obtained with the kinematic viscosity v = 10~ 4 . However, similar results are 
obtained for v — 10~ 5 . In order to establish the importance of the component of the tidal 
potential with m = 3, we consider cases for which the full tidal potential is used, only the 
m = 3 component of the tidal potential is used and finally for which the full tidal potential 
with the m = 3 component removed is used. By studying the time-evolution of the mode 
amplitudes in these cases, the importance of the m — 3 component of the tidal potential 
could be confirmed. Note that because the defined mode amplitudes involve integrals over 
the disk radius, the significance of their absolute values is unclear because of cancellation 
effects. However, they may still be used to relate to the growth rates of linear modes. 

The time-evolution of the mode amplitudes Sk,i for q = 0.1 is plotted in Fig. [2U The 
upper panel gives results for the case when only the component of the tidal potential with 
r7i = 3 is retained. The lower panel gives results corresponding to the case when the full 
tidal potential is used. In both cases the mode with k = 2 and I — 3 grows at the same 
rate as the mode with k = 1 and I — 0, the latter mode being responsible for the disk 
eccentricity. Although it is difficult to demonstrate the causal nature of this connection, 
this behaviour would certainly be expected from the discussion of mode coupling given 
above. We also find that the final state and magnitude of the eccentricity of the disk is 
similar in all cases for which the disk is unstable to becoming eccentric. However, the 
growth rate is measured to be smaller when the full potential is used. Fig. |2"21 gives similar 
plots to those shown in in Figure [5T] but for q = 0.2. Note that, as expected from the 
mode-coupling theory, there is a very close correspondence between the mode with k = 2 
and I = 3 and the mode with k = 1 and I = when only the m = 3 component of the tidal 
potential is used for this mass ratio. 
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7.4. Effectiveness of the component of the tidal potential with m = 3 

In order to determine whether the to = 3 component of the tidal potential is essential for 
producing eccentricity growth, we performed simulations using the full tidal potential with 
the to = 3 component removed. The time-evolution of the mode amplitudes Sk,i is shown 
in Fig. [221 for q = 0.1 and q = 0.2. When q = 0.1 little growth of the disk eccentricity is 
seen confirming the importance of the m = 3 component. However, when q = 0.2 some 
growth is seen but at a significantly slower rate than when the full potential is used. 

Although the results discussed above indicate the importance of the to — 3 component 
of the tidal potential, they also indicate that it is not always needed. From Fig. [23] it is 
apparent that although eccentricity growth is very small when q = 0.1, it is still present 
but reduced when q = 0.2. The time required to attain saturation is increased from about 
30 to 150 orbits in that case. 

As in the other cases the mode with k = 2 and I = 3 grows at the same rate as the 
mode with k = 1 and I = indicating that this is still responsible for causing eccentricity 
growth. Also a response with k = 3 and I = 3 is still present even though it cannot be 
driven by the to = 3 component of the potential. The most natural explanation is that this 
response is produced by a forcing originating from a coupling between the response to the 
dominant to = 2 tide and the to = 1 tide. As before that can then couple to the mode 
producing eccentricity growth. Thus in this case one has replaced the (3,3), (1,0) pairwise 
coupling by a (2, 2), (1, 1), (1, 0) triple coupling. In so doing the strength is reduced by a 
factor q making the process less significant in the smaller mass ratio case as is seen in the 
simulations. This process thus contributes more effectively for larger q. On the other hand 
the 3 : 1 resonance which enables effective angular momentum transport between disk and 
companion is pushed towards the Roche lobe requiring further disk expansion and resulting 
in the process becoming increasingly nonlinear. 

8. Comparison with linear theory 

In an attempt to understand the linear phase of the simulations, during which the eccen- 
tricity is small and grows exponentially in time, we analysed the full set of linear equations 
governing small perturbations, with azimuthal mode number m = 1, of a steady, axisym- 
metric disk. (This is in addition to the approach mentioned in section B~2l in which approx- 
imations are made that are appropriate for slowly evolving global eccentric modes.) For 
the purposes of comparison with the standard model in which a reflecting inner boundary 
condition is used, the basic state of the linear calculation is a non-accreting solution with 
zero radial velocity. 

This approach yields appropriate solutions corresponding to slowly precessing eccentric 
modes of the disk, of which the 'lowest-order' mode with the simplest radial structure 
might be expected to be the most relevant. The precession rate can be positive or negative 
depending on the mass ratio of the binary, the value of H/r, the surface density profile 
and the boundary conditions. While it is possible, using this approach, to account for a 
number of aspects of the simulation results, the sensitivity of the precession rate to the 
exact surface density profile and the boundary conditions is a matter of concern. 
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Such modes also grow or decay slowly as a result of viscosity. It is known that shear 
visco sity can render modes unstable through the mechanism of viscous overstability (jKato 
1978), in which some of the energy being drawn by viscous stresses from the differential 
rotation of the disk is diverted into growing oscillations. Viscous overstability has been 



identified as an important issue in the theory of eccentric disks (|Ogilvie 



2001). A local 



analysis in the case of an isothermal two-dimensional disk with constant kinematic shear 
viscosity and no bulk viscosity indicates overstability for radial wavelengths > 9H and a 
maximum local g rowth rate of approxim ately 0.034^/iJ 2 for a radial wavelength of ap- 
proximately 13-ff |Latter fc Qgi lvic 2006). This behaviour is confirmed in our global linear 
calculations. For sufficiently thin disks, growing modes are obtained. However, for the stan- 
dard model with H/r = 0.05, the preferred local wavelength of the viscous overstability is 
as long as 0.66r and in fact we do not find overstable global modes. In addition, we confirm 
by this method that the inclusion of a bulk viscosity with t'buik — 2 x 10 -5 makes a negative 
contribution to the growth rate that is quite negligible in comparison with those found in 
the simulations. Therefore viscous overstability is not relevant to the standard model but 
might be an important issue when thinner disks are considered. 

It is possible within this linear approach to model the nonlinear mode-coupling process 
by adding terms that correspond to a localised growth of the eccentric perturbations in 
the vicinity of the 3: 1 resonance. Rather than loc alising the growth in the form of a delta 
function as done by iGoodchild fc Qgilvid (|2006l ). this effect can be represented using a 
Gaussian or Lorentzian function of radius centred on the resonance and with a width 
appropriate to a Lindblad resonance broadened by pressure. Including these terms allows 
us to obtain growing eccentric modes with growth rates comparable to those obtained in 
the simulations. However, it is again found that a detailed comparison is difficult because 
of the sensitivity of the results to the exact surface density profile and the treatment of the 
resonance. 



According to the angular momentum equation in the case of a uniform shear viscosity, 
the surface density in the non-accreting basic state should de viate from the initia l r -1 / 2 



profile only as a result of the tidal torque applied to the disc. iPapaloizou fc Pringld (|1977l ) 
provide a method to calculate the tidal torque per unit radius, which gives a result propor- 
tional to the viscosity of the disk. When this is balanced with the viscous torque it predicts 
a surface density profile that is independent of the viscosity and depends only on the mass 
ratio q. Unfortunately this approach does not give accurate agreement with the surface den- 
sity profile shown in Fig. 01 indicating that the angular momentum balance is n onloc al and 
involves propagating density waves not considered by IPapaloizou fc Pringlc (1977) . Th e 



standard model disk is therefore smaller than predicted by IPapaloizou fc Pringld (|19770 . 
We note that a two-dimensional isothermal disk may allow radial wave propagation to a 
greater extent than more realistic models. It may also be relevant that SPH simulations are 
less likely to allow these waves to propagate and may therefore produce larger disks that 
are more susceptible to eccentric instability. Clearly the tidal truncation of disks needs to 
be better understood in order to make accurate predictions of the growth and precession 
rates of eccentric modes. 
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9. Effect of a mass-transfer stream 
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To make better contact with possibly more realistic cases with mass overflow from the 
secondary we have performed some studies where we modify the outer boundary condition 
to allow for mass inflow through the L\ point. The model parameters are dealt with in the 
same way as before. We adopt the mass ratio q = 0.1 which allows a scaling to be made 
such that binary parameters match very closely those of OY Car which has a primary of 
mass Mi — 0.685M Q and P rb = 0.063 d . The physical separation of the two stars is 
thus ctbin = O.627i?0. Mass is injected into the system at a constant rate. Note that in 
the case of locally isothermal disk models with constant kinematic viscosity as adopted 
here, the independence of the dynamics of the flow, and in particular the final values of 
eccentricity and disk radius in a quasi-steady state, to the surface density scale means that 
the results are independent of the magnitude of the inflow rate! This may thus be scaled 
to a value appropriate for cataclysmic variables, typically M = 10~ 9 AfQ/yr. However, in a 
more realistic disk model that includes suitable radiation and energy transport mechanisms 
and allows for the possibility of time-dependent outbursts, this would not be the case. 

In the simulations presented here the Roche lobe of the primary is as empty as possible 
initially, i.e. the density is constant and equals the floor value. 

The first set of simulations were performed with the standard viscosity v = 10~ 5 in 
dimensionless units which would scale to a value of 2.2 x 10 13 cm 2 /s for the OY Car system. 
The inner boundary condition is open to allow for accretion of mass onto the primary star 
and enable a steady state to be reached. We have performed simulations with both the 
'viscous outflow' and the 'free outflow' conditions as described in section [2~3l above. From 
Fig. [2U it is clear that the standard viscosity value used does not result in an eccentric 
disk. Even an increase of the viscosity by a factor of 3 does not produce an eccentric disk 
in the 'free outflow' case as the final (blue) curve in Fig. [M] indicates. In this latter case 
the disk's radius is increased to r& ~ 0.48 due to the spreading action of viscosity but the 
disk nevertheless does not switch to an eccentric mode. The apparently large initial disk 
radii in these models are explained by the fact that the disk is initially empty and the 99% 
radius is contributed to by the material that is infalling from the L\ point. 

We then ran models with a larger viscosity v = 10~ 4 using the 'viscous outflow' condi- 
tion at the inner boundary. The models were started from both an empty Roche lobe and 
a standard initial disk as in the previous simulations. In these cases the disk attained an 
eccentric mode. In both cases an identical final state is reached, however much more quickly 
for the model with an initial disk. The evolution of the total mass in the disk (within the 
primary's Roche lobe) is displayed in the top left panel of Fig. [55] The two cases show a 
complementary behaviour: in the case with an initial disk the total mass drops towards 
the final equilibrium (at around logmdisk = — 4) while for the empty case it rises initially, 
overshoots the equilibrium point and settles to the equilibrium value roughly 100 orbits 
later. Additional disk properties such as the disk radius, eccentricity and precession rate 
are all independent of the initial condition as indicated by the other panels in Fig. 1251 
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9.1. Disk dynamics 

The dynamical influence and importance of the stream inflow can now be studied (irre- 
spective of initial conditions) through a direct comparison with the corresponding model 
without a stream; the results are summarized in Table [3J The mean disk eccentricity 
and the disk radius rd in the final equilibrium phase are very similar in both cases, with an 
expected smaller disk radius and slightly smaller eccentricity due to the incoming relatively 
low angular momentum disk material. The surprising difference lies in the magnitude of the 
rate of disk precession. For the model without the stream the disk precesses in a retrograde 
sense quite rapidly, performing a full revolution in about 17 binary orbits. The impact of 
the stream is such as to slow down the precession such that it takes about 140 orbits to 
make a full revolution. 



Boundary condition 


ea 


rd 




No stream inflow 


0.56 


0.651 


-0.060 


With stream inflow 


0.55 


0.610 


-0.0072 



Table 3. Comparison of two models with and without stream inflow, all other parameters 
being identical. The viscosity is v = 10~ 4 and at the inner boundary 'viscous outflow' 
conditions are applied. The model that has no stream corresponds to the second model in 
Tabled] 

The surface density distribution of the model with stream inflow and an initial disk is 
shown during the quasi-stationary phase at 4 times separated by 10 binary orbits. From 
Fig. we can directly infer the slow retrograde precession of the disk which is also con- 
firmed independently by a phase analysis of the m = 1 Fourier mode at the specific radius 
r = 0.37. 

The change of the precession rate by about Azb = 0.52P o ~^ is also apparent from Fig. 1271 
where a quasi-stationary model with stream inflow (from Fig. I26[) is continued after the 
stream is switched off. Clearly visible is the immediate transition to the state with no 
stream. The new average precession rate corresponds closely to the case with no stream 
inflow. Quite obviously, the disk precession rate is sensitive to modifying conditions at the 
outer boundary as well as at the inner boundary as we have seen above. Additional runs 
with cooler disks (h = 0.03 and 0.04) demonstrate that the transition from prograde to 
retrograde disk precession occurs in the case with stream inflow for much larger values 
around between h — 0.04 and 0.05, see upward triangles in Fig. [TJ] On the other hand the 
disk radius and eccentricity are more robust. 

10. Conclusions 

In this paper we have performed grid-based simulations of eccentric disks in close binary 
systems. Adopting a grid-based rather than a particle-based method enables consideration 
of a wider and more appropriate range of physical parameters. We have been able to 
consider disk aspect ratios in the range 0.01 < h < 0.05 and dimensionless viscosities 
in the range 3.3 x 10~ 6 < v < 10~ 4 , whereas particle-based methods arc limited to the 
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up per limits of these ranges. We note the relationship between v and the a parameter 
of IShakura fc Sunvaevl (|1973T ) a = vjh 2 . Thus v = 10~ 4 corresponds to a ~ 0.1 for 



h = 0.03. We have also considered the effects of different inner boundary conditions, the 
effects of a mass-transfer stream and considered and demonstrated numerical convergence. 
However, due to heavy computational demands the grid-based simulations carried out 
in this first survey have been limited to two dimensions and simply adopted a locally 
isothermal equation of state. We summarize the main results of our simulations below. 



10.1. Instability to the formation of an eccentric disk 

We find instability to the formation of an eccentric disk in the mass range 0.05 < q < 0.3. 
Explorative calculations have shown instability even for larger mass ratios up to q = 1 in 
these two-dimensional simulations, see Fig. The most unstable cases are those with a 
reflecting inner boundary, large viscosity, small aspect ratio and no mass-transfer stream. 
An initially non-eccentric disk quickly enters a linear phase of exponential growth of an 
unstable mode. 

For q — 0.1 we found growth times ranging between 162 and 15 binary orbits for 
v in the range 3.3 x 10~ 6 — 10~ 4 with h in the range 0.02 — 0.05. A quasi-stationary 
precessing eccentric disk is attained in all cases in about four growth times. In the most 
unstable cases (large v and small h) this corresponds to about 60 binary orbits or about 6 
days for superoutbursting cataclysmic binaries. This is consistent with reported superhump 
development times of a few days. 

The mean eccentricities of the quasi-stationary disks are typically in the range 0.3 — 0.6 
when there is no mass-transfer stream and an inn er reflecting boundary . This is in the 



up per part of the ran ge 0.38 ±0.1 found for OY Car lHessman et alj (| 1992T ) and for IY Ma 

by 



Rolfe et al 



(|200lD . 

However, it is important to note that results are sensitive to the inclusion of a mass- 
transfer stream and the treatment of the inner boundary condition. For the largest viscosi- 
ties, use of an inner outflow condition rather than a reflecting boundary condition reduces 
the final mean eccentricity from 0.5 to 0.3. Incorporation of a mass-transfer stream does not 
affect the growth of eccentricity for v = 10~ 4 but we are able to see a transition to stability 
to formation of an eccentric disk at about v = 3 x 10~ 5 in this case independently of the in- 
ner boundary condition. This is also consistent with the idea that the effective viscosity and 
consequently the mass-transfer rate through the disk is larger during the outburst ph ases 
of cataclysmic binaries during which supcrhumps are reported ([Patterson et al.l l2005). 



10.2. Quasi-steady state of the eccentric disk 

Once the simulated disks go through an instability that results in them becoming eccentric, 
they typically attain a characteristic quasi-steady state. This was illustrated in Fig. [TU1 
Viewed from an inertial frame the secondary has a close approach to the slowly precessing 
disk apocentre approximately once every orbit, pulling out a tidal tail. This later impacts 
back on the disk which is not strongly affected by the secondary until the next close 
approach. Thus almost the entire tidal dissipation may occur as a result of these events. In 
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this context we note that tidal dissipation has been proposed as the source of the luminosity 
variations during the main superhump phase while it has been proposed that effects due 
to mo dulation by the mass-transfer stream become significant durin g the late superhump 



phase ([Hessman et al 



f992 



Rolfe et al 



2001 



Patterson et al 



2005). 



10.3. Instability mechanism 

We have investigated the mechanism responsible for the excitation of the disk ecce ntricity 
and fou nd this to be basically consistent with the mode-coupling mechanism of 



Lubow 



(|1991af) . Neglecting the slow precession of the disk, viewed from the inertial frame the 
eccentricity of the disk may be considered to arise from a stationary surface density per- 
turbation with azimuthal mode number m = 1. In addition the disk responds directly to 
the tidal potential of the companion producing a surface density response that corotates 
with it and consists of a sum of Fourier components with m — 1,2,3,... with amplitudes 
that scale with the mass ratio, q. Nonlinear coupling between the response to the tidal 
potential and the m = 1 mode causes the excitation of secondary density waves with pat- 
tern speeds muj/{m — 1), m = 2,3, 4, ... which exceed that of the binary and are thus 
associated with a tidal energy and angular momentum transfer to it. The back reaction 
on the m = 1 mode causes the eccentricity to grow. To work effectively the secondary 
waves should be associated with an inner Lindblad resonance. The most favourable case 
corresponds to m = 3 and the disk gas orbiting with Q = 3a;, at a 3:1 commensurability. 
In this case the Fourier component of the tidal potential with m — 3 plays a key role. 

The results of simulations with q — 0.1 supported this view with growth being almost 
absent when the m = 3 component was removed. However, for q = 0.2 growth still occurred 
in this situation albeit at a reduced rate. This is apparently because higher-order couplings 
generated through the other Fourier components result in an effective forcing with m = 3. 
Being of higher order, this is more effective at larger q. But it should be noted that as 
q increases, the 3 : 1 commensurability is driven towards the Roche lobe thus becoming 
increasingly nonlinear and affected by factors such as the outer boundary condition and any 
stream inflow. As we indicate below, it is possible that dissipative effects not included in our 
two-dimensional simulations play a more important role. For these reasons we concentrated 
on the case q = 0.1 in this article. 



10.4. Disk precession rates and the dependence on disk aspect ratio 

The main determinants of the disk precession rate were found to be the disk aspect ratio 
and the impact of a mass-transfer stream if present. Results for 0.01 < h < 0.05 for q = 0.1 
are summarized in Fig. 1141 

The quasi-steady disk precession rate was found to decrease from positive values 
at small h (prograde precession) to attain negative values at larger h. This is because, as 
expected, the pressure contribution tends to produce retrograde precession. When no mass- 
transfer stream was present the transition from prograde to retrograde occurred for h ~ 
0.025. We also determined the precession rates by independently solving the linear mode 



Kley at al.: Eccentric disks in binary stars 



35 



problem using the average surface density profile and in view of the strong nonlinearity 
present obtained surprising agreement with h at the transition being ~ 0.03. 

An approximate parabolic fit to the data points for the case with no mass-transfer 
stream was f ound to be w& = 0.022 — 39h 2 . This indicates that roj ~ 0.02 as indicated by 



observations (Patterson et al 



2005) only for h ^ 0.01. In view of the fac t that values for h 



estimated from accretion d isk modelling fall in the range 0. 01 < h < 0.03 ( Smith et al 



Goodchild fc Q gilvie 2006), this is rather small. Note that 



Goodchild fc Qgilviel (|2006l) also 



2007 



obtained similarly small estimates from their one-dimensional linear formalism. 

However, the simulations for which the disk was impacted by a mass-transfer stream 
show the transition for h w 0.04 and indicate the possibility of obtaining the observed 
precession rates for values of h estimated from accretion disk modelling. Thus our results 
demonstrate the importance of the mass-transfer stream in this context as well as in the 
context of obtaining a transition from instability to stability at a reasonable value for the 
disk viscosity. On the other hand, it might be argued that the significance of the mass- 
transfer stream would be reduced during a superout burst, when the accretion rate through 
the disk exceeds the mass-transfer rate from the companion. 

In addition we comment that there is sensitivity to the inner boundary condition such 
that changing from a reflecting condition to a dissipative or free outflow condition also acts 
to make the precession frequency move in the prograde direction. 



10.5. Comparison with particle-based simulations 



Smith et al 



(2007). 



The most recent relevant SPH simulations have been performed by 
Both two- and three-dimensional simulations have been considered. Several important is- 
sues should be noted. Such simulations are inevitably highly diffusive (a > 0.1) and thick 
(h > 0.05) so that parameter studies at smaller h and smaller a of the type considered 
here cannot be carried out. Comparison has to be limited to cases with large v and h. The 
SPH simulations are performed with ~ 10 5 particles and also with a mass-transfer stream. 
For a disk with h = 0.05 one can estimate that there is only about one smoothing length 
per scaleheight in a three-dimensional simulation so that the vertical thickness is poorly 
resolved. The three-dimensional SPH simulations tend to produce a weaker phenomenon 
than that reported here with almost no activity for q > 0.2 and mean disk eccentricities 
at the small end of the range implied by observation. In addition very long growth times 
approaching thousands of orbits are generally reported. 

On the other hand 



Patterson et al 



( 20051) indicate that the super hump phenomenon 



occurs for q < 0.35 and that it is associated with high mass-transfer events during super- 
outbursts and a short development time requiring growth times ~ 20 binary orbits (see 
above). It may also occur in a quasi-steady manner in systems such as nova-like variables 
which are also believed to have a high mass-transfer rate which in this case is steady. This 
is also consistent with our finding that stability to eccentricity generation occurs at low 
viscosities and accordingly low mass-transfer rates through the disk when a mass-transfer 
stream is present. 
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The reported two-dimensional SPH simulations are more in line with those reported 
here. This suggests that the three-dimensional case is more dissipative though the adequacy 
of the resolution in the vertical direction should also be an issue. It is likely that inclusion 
of effects arising from the vertical structure of the disk could result in processes leading 



to additional dissipation (|Papaloizoul l2005f ) but these would need to be simulated with 
adequate resolution. The two-dimensional grid-based simulations reported here are possibly 
too unstable to forming eccentric disks, particularly at larger q. However, in those cases, the 
3:1 resonance driving the instability is pushed towards the Roche lobe where higher-order 
mode couplings and increasing nonlinearity occur and numerical studies should include a 
treatment of the vertical direction if this region is to be correctly dealt with. 

Therefore an important direction for future work will be to develop three-dimensional 
grid-based simulations which resolve the disk vertical structure and address the issue of 
the importance of possible additional dissipation modes. 
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Fig. 21. Time-evolution of the mode amplitudes Sk,i (see Eq. IT3| for q = 0.1 with a 
kinematic viscosity v = 10 plotted on a logarithmic scale. The upper panel gives results 
for the case when only the component of the tidal potential with m = 3 is retained. The 
lower panel gives results corresponding to the case when the full tidal potential is used. In 
both cases the mode with k — 2 and I = 3 grows at the same rate as the mode with k = 1 
and I = 0, the latter mode being responsible for the disk eccentricity. The growth is slower 
when the full potential is used. 
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Fig. 22. As in Figure I2T1 but for q — 0.2. Note that there is a very close correspondence 
between 52,3 and Si t o when only the m = 3 component of the tidal potential is used. 
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Fig. 23. Time-evolution of the mode amplitudes Sk,i, plotted on a logarithmic scale, for 
q = 0.1 (upper panel) and q = 0.2 (lower panel). For these cases, the tidal potential with 
the m = 3 component removed is used and the kinematic viscosity v = 10~ 4 . When q = 0.1 
little growth of the disk eccentricity or Si t o is seen. When q = 0.2 some growth occurs but 
at a significantly slower rate than when the full potential is used (see Fig. [22]) . 
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Fig. 24. Mean disk eccentricity ed and disk radius r<j as functions of time for models that 
include mass inflow through the L\ point and different inner boundary conditions and 
viscosities. For the first two plots (green and red curves) the standard viscosity 10 -5 has 
been used, while the third plot (blue curve) is for a model that has been continued from 
the first at time t 230 but with v = 3x 10~ 5 . 
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Fig. 25. Disk mass md, radius rd and eccentricity ed as functions of time for models 
including mass inflow through the L\ point with different initial conditions. The first (red) 
curve corresponds to an initially empty Roche lobe of the primary, while in the second 
(green) a standard initialized disk is used. The 'viscous outflow' condition and a viscosity 
v = f0 -4 (ten times standard) are used. 
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Fig. 26. Greyscale plot of the surface density of the disk for the model with stream inflow 
that started with an initial disk at 210, 220, 230 and 240 orbital periods. The solid (red) 
curve indicates the Roche lobe of the primary star (central red dot). The 'viscous outflow' 
condition is used at the inner boundary. 
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Fig. 27. Longitude of the disk pericentre (in the inertial frame) as a function of time. 
The solid (red) curve corresponds to a disk in a quasi-steady state with stream inflow (as 
displayed in the previous Fig. [26)) . The light dashed (green) line corresponds to a model 
where the stream is switched off at t = 240 and the simulation is continued. 



